#include <math.h>
#include <vector>
using std::vector;
	static const double fPDF[] = {
		2.3283064365386963e-10,7.450580596923828e-9,
		1.1548399925231934e-7,1.1548399925231934e-6,
		8.372589945793152e-6,0.00004688650369644165,
		0.00021098926663398743,0.000783674418926239,
		0.002448982559144497,0.0065306201577186584,
		0.015020426362752914,0.03004085272550583,
		0.0525714922696352,0.08087921887636185,
		0.10976465418934822,0.13171758502721786,
		0.13994993409141898,0.13171758502721786,
		0.10976465418934822,0.08087921887636185,
		0.0525714922696352,0.03004085272550583,
		0.015020426362752914,0.0065306201577186584,
		0.002448982559144497,0.000783674418926239,
		0.00021098926663398743,0.00004688650369644165,
		8.372589945793152e-6,1.1548399925231934e-6,
		1.1548399925231934e-7,7.450580596923828e-9,
		2.3283064365386963e-10};

class BinomialStat {
	double m_Total;
	double m_SqTotal;
	size_t m_n;
	double fRmsMax;
	vector<double> fCounts;
	vector<double> fTotal;
	vector<double> fSqrTotal;
	double fProb[33][33];
 public:
	BinomialStat() : m_Total( 0.0 ), m_SqTotal( 0.0 ), m_n( 0 ),
		fCounts(33, 0.0), fTotal(32, 0.0), fSqrTotal(32, 0.0)
		{
			memset(fProb, 0, sizeof(double)*33*33);
		}

	double rms(int i)
		{
			fRmsMax = -1.0;
			double sum = 0.0;
			double tmp = 0.0;
			double count = 0;
			for (int j = 12; j <= 20; ++j) {
				tmp = fPDF[j] - (fProb[i][j] / fCounts[i]);
				tmp *= tmp;
				sum += tmp;
				if (tmp > fRmsMax) fRmsMax = tmp;
				count++;
			}
			return ::sqrt(sum / count);
		}

	double rms()
		{
			fRmsMax = -1.0;
			double sum = 0.0;
			double tmp = 0.0;
			double count = 0;
			for (int i = 0; i < 32; ++i) {
				for (int j = 12; j <= 20; ++j) {
					tmp = fPDF[j] - (fProb[i][j] / fCounts[i]);
					tmp *= tmp;
					sum += tmp;
					if (tmp > fRmsMax) fRmsMax = tmp;
					count++;
				}
			}
			return ::sqrt(sum / count);
		}

	double rms_max() const
		{
			return ::sqrt(fRmsMax);
		}

	void insert(int i, int j)
		{
			fProb[i][j] += 1.0;
			fCounts[i] += 1.0;
			fTotal[i] += j;
			fSqrTotal[i] += j*j;

			m_Total += j;
			m_SqTotal += j*j;
			++m_n;
		}

	void clear()
		{
			m_Total = 0.0;
			m_SqTotal = 0.0;
			m_n = 0;
			fRmsMax= -1.0;

			fCounts.clear(); fCounts.insert(fCounts.begin(), 33, 0.0);
			fTotal.clear(); fTotal.insert(fTotal.begin(), 32, 0.0);
			fSqrTotal.clear(); fSqrTotal.insert(fSqrTotal.begin(), 32, 0.0);

			memset(fProb, 0, sizeof(double)*33*33);
		}

	double mean(int i) const
		{
			return fTotal[i] / fCounts[i];
		}

	double mean() const // m_n must be non-zero
		{
			return m_Total / m_n;
		}

	double variance(int i) const
		{
			const double tmp = fTotal[i];
			const double count = fCounts[i];
			return (fSqrTotal[i] - tmp*tmp / count) / count;
		}

	double variance() const
		{
			return (m_SqTotal - m_Total * m_Total / m_n) / m_n;
		}

	double sd() const
		{
			return ::sqrt( variance() );
		}

	double sd(int i) const
		{
			return ::sqrt(variance(i));
		}
};


